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ABSTRACT 

Many  background  error  correlation  (BEC)  models  in  data  assimilation  are  formulated  in  terms  of  a 
smoothing  operator  B,  which  simulates  the  action  of  the  correlation  matrix  on  a  state  vector  normalized  by 
respective  BE  variances.  Under  such  formulation,  B  has  to  have  a  unit  diagonal  and  requires  appropriate 
renormalization  by  rescaling.  The  exact  computation  of  the  rescaling  factors  (diagonal  elements  of  B)  is 
a  computationally  expensive  procedure,  which  needs  an  efficient  numerical  approximation. 

In  this  study  approximate  renormalization  techniques  based  on  the  Monte  Carlo  (MC)  and  Hadamard 
matrix  (HM)  methods  and  on  the  analytic  approximations  derived  under  the  assumption  of  the  local  ho¬ 
mogeneity  (LHA)  of  B  are  compared  using  realistic  BEC  models  designed  for  oceanographic  applications.  It 
is  shown  that  although  the  accuracy  of  the  MC  and  HM  methods  can  be  improved  by  additional  smoothing, 
their  computational  cost  remains  significantly  higher  than  the  LHA  method,  which  is  shown  to  be  effective 
even  in  the  zeroth-order  approximation.  The  next  approximation  improves  the  accuracy  1.5-2  times  at 
a  moderate  increase  of  CPU  time.  A  heuristic  relationship  for  the  smoothing  scale  in  two  and  three  di¬ 
mensions  is  proposed  for  the  first-order  LHA  approximation. 


1.  Introduction 

Modeling  of  the  background  error  correlation  (BEC) 
by  smoothing  operators  in  variational  data  assimilation 
has  recently  gained  considerable  attention,  primarily  due 
to  computational  efficiency  of  their  implementation  and 
versatility  in  approximating  anisotropic  and  inhomoge¬ 
neous  BECs  (e.g.,  Xu  2005;  Pannekoucke  and  Massart 
2008;  Mirouze  and  Weaver  2010).  In  the  framework  of 
this  approach,  the  correlation  matrix  C  of  the  BE  field 
is  approximated  by  a  positive  function  of  the  diffusion 
operator: 

D  -  Vi/V,  (1) 

where  v  is  the  diffusion  tensor,  whose  components  may 
depend  on  the  spatial  coordinates  x.  Attractive  features 
of  such  models  are  their  low  computational  cost  and 
direct  control  of  inhomogeneity  and  anisotropy  by  v 
under  the  positive-definiteness  constraint.  Among  the 
most  popular  functions  of  D  used  in  practice,  are  the 
exponential  (yielding  the  Gaussian-shaped  correlations) 
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and  the  inverse  of  a  positive  polynomial.  Numerically, 
these  functions  are  implemented  by  integrating  the  finite- 
difference  diffusion  equation  on  the  model  grid  and  have 
been  used  in  many  applications  (e.g.,  Derber  and  Rosati 
1989;  Weaver  et  al.  2003;  Di. Lorenzo  et  al.  2007;  Liu 
et  al.  2009). 

An  important  numerical  issue  arising  with  this  type  of 
covariance  modeling  is  the  necessity  to  find  a  rescaling 
transformation  R  such  that  the  diagonal  elements  of  the 
covariance  matrix  are  equal  to  the  background  error 
(BE)  variances  derived  from  the  history  of  data  assimi¬ 
lation  into  the  numerical  model.  The  direct  way  of  finding 
R  is  to  compute  the  inverse  square  root  of  the  diagonal 
elements  of  C.  The  latter  can  be  found  by  convolving  C 
with  5  functions  in  every  point  x  of  the  model  grid  and 
taking  the  results  of  these  convolutions  at  the  same  points. 
This  procedure,  however,  is  practically  unfeasible  as  the 
model  grid  size  (the  number  of  convolutions  to  be  exe¬ 
cuted)  may  often  exceed  M  ~  106-107. 

In  this  study  we  compare  the  numerical  efficiency  of 
several  methods  of  estimating  the  diagonal  elements  of  a 
symmetric  positive-definite  matrix  C  defined  by  its  ac¬ 
tion  on  an  arbitrary  vector  and  propose  an  algorithm  for 
estimating  R  under  the  locally  homogeneous  approxi¬ 
mation  (i.e.,  slowly  varying  smoothing  scales  of  C).  The 
manuscript  is  organized  as  follows.  In  section  2,  we  briefly 
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overview  the  diagonal  estimation  methods  with  an  em¬ 
phasis  on  the  possibility  to  exploit  the  analytic  formulas, 
available  for  the  considered  functional  forms  under  the 
assumption  of  spatial  homogeneity.  In  particular,  an 
extension  of  the  method  of  Purser  et  al.  (2003)  for  esti¬ 
mating  the  Gaussian  kernel  diagonal  is  proposed  for 
higher  dimensions.  In  section  3,  we  describe  the  setup  of 
numerical  experiments  with  anisotropic  inhomogeneous 
correlation  model  and  document  the  results.  Summary 
and  conclusions  complete  the  paper. 

2.  Renormalization  methods 

a.  Monte  Carlo  technique 

This  method  originates  from  a  large  family  of  stochastic 
algorithms  used  for  estimating  elements  and  traces  of 
extra-large  matrices  (e.g.,  Girard  1989;  Hutchison  1989; 
Dong  and  Liu  1994).  Weaver  and  Courtier  (2001)  were 
among  the  first  to  use  this  approach  in  geophysical  ap¬ 
plications  for  estimating  the  diagonal  of  the  Gaussian- 
shaped  BEC  operators. 

The  underlying  idea  is  to  define  an  ensemble  of  K  ran¬ 
dom  vectors  s*  on  the  model  grid  and  perform  compo¬ 
nentwise  averaging  of  the  products  s  —  Cs  according  to 
the  following  formula: 

d(x)  =  s  G  s0s  O  s,  (2) 

where  the  overbar  denotes  averaging  over  the  ensemble, 
and  Oand  0  stand  for  the  componentwise  multiplication 
and  division  of  the  vectors,  respectively.  Simple  consid¬ 
erations  show  that  when  all  the  components  of  s  have 
identical  5-correlated  distributions  with  zero  mean,  con¬ 
tributions  to  d  from  the  off-diagonal  elements  tend  to 
cancel  out,  and  d  converges  to  d  =  diag(C)  as  K  oo. 
More  accurately,  the  squared  relative  approximation  error: 

e2(x)  =  (d  -  d)2/d2  (3) 


way  to  construct  such  an  ensemble  is  to  draw  the  vectors 
s*  from  the  columns  of  the  M  X  M  Hadamard  matrix 
(HM),  which  span  the  model’s  state  space  (see  appendix 
A  for  more  details  on  the  HM  theory). 

Although  it  is  not  proven  yet  that  Hadamard  matrices 
can  be  constructed  for  an  arbitrary  M,  very  efficient  re¬ 
cursive  algorithms  for  generating  HM  columns  do  exist 
for  Ms ,  whose  factorization  involves  only  prime  num¬ 
bers  not  exceeding  100.  Since  the  exact  convergence 
(which  can  be  achieved  at  M  iterations)  is  never  needed 
in  practice,  it  is  not  even  necessary  to  draw  s*  from  the 
Hadamard  matrix,  whose  dimension  exactly  coincides 
with  the  state  space  dimension:  if  generation  of  s*  is  im¬ 
possible  because  of  some  odd  value  of  M  (e.g.,  1004),  it 
can  always  be  replaced  by  a  slightly  larger  number  (e.g., 
1008),  and  sk  can  be  defined  as  the  first  1004  components 
of  the  1008-dimensional  Hadamard  vectors.  In  the  nu¬ 
merical  experiments  described  in  the  next  section  it  is 
shown  that  such  modification  does  not  affect  the  method’s 
convergence  during  the  first  several  hundred  iterations. 

c.  Locally  homogeneous  approximations 

An  alternative  approach  to  the  diagonal  estimation 
employs  a  priori  information  on  the  structure  of  C.  Con¬ 
sider  a  homogeneous  {v  =  const)  case  and  assume  that 
the  coordinate  axes  are  aligned  along  the  eigenvectors 
of  the  diffusion  tensor,  whose  (positive)  eigenvalues  are 
A?,i  =  1, ... ,  n  and  n  is  the  number  of  dimensions  of  the 
physical  space  (1,  2,  or  3).  Then  the  matrix  elements  of 
the  two  types  of  BEC  operators  that  are  considered  here 
can  be  written  down  explicitly  (see  appendix  B)  as 


C,(*,y) 


exp(D/2)  =  dexp 


(4) 


C2(x,y) 


(I  -  D/2rn)_m 


2*-'r(Sy 


(5) 


where 


is  inversely  proportional  to  the  ensemble  size  K.  In  other 
words,  one  may  expect  to  achieve  10%  accuracy  at  the 
expense  of  approximately  100  multiplications  by  C  if  the 
first  ensemble  member  gives  a  100%  error.  This  estimate 
may  seem  acceptable  since  in  geophysical  applications 
the  BE  variances  are  usually  known  with  limited  preci¬ 
sion  and  approximating  the  diagonal  with  5%-10%  er¬ 
ror  seems  satisfactory. 

b.  The  Hadamard  matrix  method 

The  Monte  Carlo  (MC)  technique  was  developed  fur¬ 
ther  by  Bekas  ct  al.  (2007),  who  noticed  that  the  method 
may  converge  to  d  in  the  finite  number  of  iterations  M  if 
the  ensemble  vectors  are  mutually  orthogonal.  An  easy 


p  =  v(x-y)T"-,(x  ~  y) 

is  the  distance  between  the  correlated  points,  measured 
in  terms  of  the  smoothing  scales  A,-;  d  =  (2n-)“rt/2fl'1  are 
the  (constant)  diagonal  elements  of  C^;  Ct  =  IIAi  = 
\/deti/  is  the  diffusion  volume  element;  I  is  the  identity 
operator;  K  and  V  denote  the  Bessel  function  of  the  sec¬ 
ond  kind  and  the  gamma  function,  respectively;  m  is 
a  positive  integer  (an  approximation  order  of  by  C2); 
s  —  m  —  nl 2;  and  p  =  \/2mp.  The  parameter  m  can  also  be 
interpreted  as  the  number  of  “time  steps”  used  in  discrete- 
time  integration  of  the  corresponding  diffusion  equation 
by  the  implicit  scheme  (see  appendix  B). 
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When  v  varies  in  space,  (4)-(5)  are  no  longer  valid, 
and  the  diagonal  elements  d  depend  on  x  and  on  the  type 
of  the  BEC  operator.  However,  if  we  assume  that  v  is 
locally  homogeneous  (LH)  (i.e.,  varies  in  space  on  a  typical 
scale  L ,  which  is  much  larger  than  A;)  the  diagonal  ele¬ 
ments  d(x)  can  be  expanded  in  the  powers  of  the  small 
parameter  e  =  A/L,  where  A  is  the  mean  eigenvalue  of 
v/F.  The  zeroth-order  LH  approximation  term  (LHO)  is 
apparently 

d°(x)  =  (2ir)_n/2fl(x)_1  (6) 


exp(A/2)~l  +  1a,  (9) 

and  represent  the  second  term  in  the  round  brackets  of 
(7)  as  follows: 

V  •  divh  =  -A  trh  +  V  •  divh',  (10) 

n 

where  h'  is  the  traceless  part  of  h.  On  the  other  hand,  if 
we  neglect  the  divergence  of  h',  (7)  can  be  rewritten  in 
the  following  form: 


because  for  infinitely  slow  variations  of  v  {L  ®)  the 
normalization  factors  must  converge  to  the  above  ex¬ 
pression  for  the  constant  diagonal  elements  d .  It  is  note¬ 
worthy  that  (6)  was  found  to  be  useful  even  in  the  case  of 
strong  inhomogeneity  e  >  1.  In  particular,  Mirouze  and 
Weaver  (2010)  found  that  such  an  approximation  pro¬ 
vided  10%  errors  in  a  simplified  ID  case. 

The  accuracy  of  (6)  can  formally  be  increased  by  con¬ 
sidering  the  next  term  in  the  expansion  of  the  diagonal 
elements  of  Cu<  The  technique  of  such  asymptotics  has 
been  well  developed  for  the  diagonal  of  the  Gaussian 
kernel  (4)  in  Riemannian  spaces  (e.g.,  Gusynin  and 
Kushnir  1991;  Avramidi  1999).  More  recently,  the  ap¬ 
proach  was  considered  by  Purser  (2008a,b)  in  the  at¬ 
mospheric  data  assimilation  context.  Application  of 
this  technique  to  the  diffusion  operator  (1)  in  flat  space 
yields  the  following  asymptotic  expression  for  the  diag¬ 
onal  elements  of  C1  in  the  local  coordinate  system,  where 
v(x)  is  equal  to  the  identity  matrix  and  D  takes  the  form  of 
the  Laplacian  operator  A: 

-  5trh  -  n(it,h  +  vdi,h)] 

+  0(<?).  (7) 

Here  h  is  a  small  (|h|  ~  e)  correction  to  *rin  the  vicinity  of 
x.  Note  that  terms  in  the  round  brackets  have  the  order 
0(c3)  because  each  spatial  differentiation  adds  an  extra 
power  of  e 

The  asymptotic  estimate  (7)  involves  second  derivatives 
that  tend  to  amplify  errors  in  practical  applications  when  e 
may  not  be  small.  Therefore,  using  (7)  in  its  original  form 
could  be  inaccurate  even  at  a  moderately  small  value  of  e. 
To  increase  the  computational  efficiency,  it  is  also  desir¬ 
able  to  formulate  the  first-order  approximation  as  a  linear 
operator,  acting  on  d°(x).  Keeping  in  mind  that  |h|  ~  e,  we 
can  utilize  the  following  relationships: 

d°(x)  =  (2i7)“n/2a(x)“1  ~ (2ir)~"/2 -  itrh) ,  (8) 


C,(x,x) 


where 


yn 


1  l_ 
6  +  3  n 


(12) 


Taking  (8)-(9)  into  account  and  replacing  A  by  D,  we 
finally  get  the  desired  ansatz  for  the  first-order  approx¬ 
imation  (LH1)  of  the  diagonal: 

d1  =  exp(vnf)d°-  (13) 


The  relationship  (13)  was  derived  by  Purser  et  al.  (2003) 
for  one-dimensional  case  (yi  =  0.5)  and  tested  by  Mirouze 
and  Weaver  (2010),  who  reported  a  significant  (2-4  times) 
improvement  of  the  accuracy  in  ID  simulations. 

It  is  likely  that  the  estimate  similar  to  (13)  can  also  be 
obtained  for  C2,  possibly  with  a  different  coefficient  yn. 
We  assume,  however,  that  yn  may  not  differ  too  much 
from  yn  given  similarity  in  the  shapes  of  the  correlation 
functions  (4)-(5)  (Fig.  1).  Furthermore,  because  of  the 
approximate  nature  of  (13),  the  best  representation  of 
d(x)  in  realistic  applications  may  be  achieved  with  a 
value  of  yn  significantly  different  from  the  one  given  by 
(12).  For  that  reason,  in  the  numerical  experiments  we 
adopt  a  more  general  form  of  (13),  assuming 


d}(x) ~ exp[yD/2]d?(x);  d> (x) ~[l- yD/4]-2d°(x), 

(14) 


and  investigate  the  dependence  of  the  respective  approx¬ 
imation  errors  {e\j}  on  the  free  parameter  y. 


3.  Numerical  results 

To  assess  the  efficiency  of  the  methods  outlined  in  the 
previous  section,  two  series  of  numerical  experiments 
with  realistically  inhomogeneous  BEC  models  are  per¬ 
formed.  In  the  first  series  we  test  the  methods  in  the  2D 
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FlO.  1.  Correlation  functions  corresponding  to  the  homogeneous 
operators  (4)-(5)  with  identical  decorrelation  scales. 


case  with  the  state  vector  having  the  dimension  of  sev¬ 
eral  thousand.  In  the  second  series,  the  LHO  and  LH1 
techniques  are  examined  in  a  realistic  3D  setting  with  a 
state  space  dimension  of  — 106. 

a.  Experimental  setting  in  2D 

The  state  space  is  described  by  scalar  functions  defined 
on  the  orthogonal  curvilinear  grid  of  the  Navy  Coastal 
Ocean  Model  (NCOM;  Martin  et  al.  2009)  set  up  in  the 
Monterey  Bay  (Fig.  2).  The  number  M  of  grid  points 
(dimension  of  the  state  space)  was  3438.  A  vector  field 
v(x)  was  used  to  generate  the  diffusion  tensor  as  follows. 
The  smaller  principal  axis  A2  of  y/v  is  set  to  be  orthog¬ 
onal  to  v  with  the  corresponding  “background”  length 
scale  A 2  =  35,  where  8(x)  is  the  spatially  varying  grid 
step.  The  length  of  the  larger  axis  A  i  is  set  to  be  equal  to 
max(l,  y/\y\/v)\v  where  v  is  a  prescribed  threshold  value 
of  |v|.  If  v  is  a  velocity  field,  then  a  structure  like  this 
simulates  enhanced  diffusive  transport  of  model  errors  in 
the  regions  of  strong  currents  on  the  background  of  iso¬ 
tropic  error  diffusion  with  the  decorrelation  scale  A2. 

In  the  2D  experiments,  the  vector  field  v  is  generated 
by  treating  bottom  topography  fc(x)  (Fig.  2)  as  a  stream- 
function.  The  threshold  value  v  was  taken  to  be  one-fifth 
of  the  rms  variation  of  |V/j|  over  the  domain. 

All  the  experiments  described  in  sections  3b-3e,  are 
performed  using  the  BEC  models  (4)-(5)  with  the  para¬ 
meters  n  =  m  =  2.  Composite  maps  of  five  columns  of  the 
corresponding  BEC  operators  are  shown  in  Figs.  2a, b. 
The  diffusion  operator  (1)  is  constrained  to  have  a  zero 
normal  derivative  at  the  open  and  rigid  boundaries  of 
the  domain  in  both  2D  and  3D  experiments. 

Numerically,  the  action  of  the  Gaussian-shaped  BEC 
operator  Cj  on  a  state  vector  y0  was  evaluated  by  ex¬ 
plicitly  integrating  the  corresponding  diffusion  equation 
y,  =  D/2y  for  the  virtual  “time  period”  defined  by  v 


starting  from  the  “initial  condition”  y0.  The  minimum 
number  of  “time  steps”  required  for  the  scheme’s  stability 
in  such  setting  was  5256.  The  action  of  C2  was  computed 
by  solving  the  system  of  equations  (I  -  D/4)2y  =  y0  with 
a  conjugate  gradient  method.  The  number  of  iterations, 
required  for  obtaining  a  solution,  varied  within  2000- 
2500.  To  make  the  shapes  of  the  C1  and  C2  compatible 
(Fig.  1),  the  diffusion  tensor  in  C2  was  multiplied  by  8/77 
(Yaremchuk  and  Smith  2011). 

The  exact  values  d(x)  of  the  diagonal  elements  are 
shown  in  Fig.  2  (right  panel).  Their  magnitude  appears 
to  be  lower  in  the  regions  of  “strong  currents”  (large  v), 
as  the  corresponding  S  functions  are  dispersed  over  larger 
areas  by  diffusion.  The  d(x)  are  higher  near  the  bound¬ 
aries  because  part  of  the  domain  available  for  dispersion 
is  screened  by  the  condition  prescribing  zero  flux  across 
either  the  open  or  rigid  boundary. 

b.  Monte  Carlo  technique 

The  MC  method  is  implemented  in  two  ways:  in  the 
first  series  of  experiments,  the  components  of  s*  are  taken 
to  be  either  1  or  -1  with  equal  probability;  in  the  second 
series  they  are  drawn  from  the  white  noise  on  the  interval 
[-1,  1],  The  residual  error  e  is  computed  using  (3).  In 
both  series  the  rates  of  reduction  of  e  with  iteration  k 
were  similar  and  closely  followed  the  \fk  law  (thin  solid 
line  in  Fig,  3c). 

Figure  3a  shows  the  distribution  of  e(x)  after  60  iter¬ 
ations  of  the  MC  method  with  the  C2  BEC  model.  The 
estimate  is  apparently  affected  by  sampling  noise,  which 
can  be  identified  by  fine  structures  at  scales  close  to  the 
grid  spacing. 

To  improve  the  accuracy,  the  MC  estimates  are  low- 
pass  filtered  with  the  corresponding  BEC  operators  at 
every  iteration.  To  optimize  the  filter,  the  diffusion  oper¬ 
ator  in  C  is  multiplied  by  the  tunable  parameter  y,  which 
effectively  reduced  the  mean  decorrelation  (smoothing) 
scale  k  =  y~m  times.  Figure  3b  demonstrates  the  result 
of  such  optimal  smoothing  (*opt  =  2.5),  which  enabled 
almost  fourfold  reduction  of  the  domain-averaged  error 
(e)  to  0.15.  The  right  panel  in  Fig.  3  shows  the  evolution 
of  (e)  with  iterations  for  the  C2  operator.  The  optimal 
smoothing  scale  reduction  factor  K~i  appears  to  follow 
the  law  k_1(/c)  ~  k~1/3  (dashed  line  in  Fig.  3c).  It  is  re¬ 
markable  that  significant  error  reduction  occurs  even  if 
the  smoothing  scale  falls  below  the  grid  spacing  k~1  <  V3. 

c.  The  Hadamard  matrix  method 

Experiments  with  the  Hadamard  scheme  are  also  done 
in  two  series.  Since  the  value  of  M  =  3438  is  not  divisible 
by  4,  it  is  hard  to  find  the  HM  of  that  dimension.  In¬ 
stead,  in  the  first  series  of  experiments,  s*  were  speci¬ 
fied  as  the  first  3438  numbers  taken  from  the  columns  of 
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the  3456-dimensional  HM,  which  can  be  easily  constructed 
from  the  12-dimensional  HMs. 

To  check  the  impact  of  ‘‘nonexact’1  dimension  on  the 
convergence,  we  artificially  increased  M  to  3456  by  re¬ 
moving  18  land  points  in  the  domain.  Experiments  with 
the  HM  of  exact  dimension  show  that  differences  in 
convergence  between  the  nonexact  and  exact  experiments 
start  to  be  visible  only  after  several  hundred  iterations. 
After  M  iterations  the  error  of  the  exact-HM  method  is 
reduced  to  the  machine  precision,  while  the  error  of  the 
first  series  of  experiments  stumbled  at  approximately  10-3 
after  1500  iterations.  This  is  consistent  with  the  18/3438  ~ 
0.5%  degree  of  nonorthogonality  of  s*.  drawn  from  the 
nonexact  HM. 

Similar  to  the  MC  method,  the  accuracy  of  HM  esti¬ 
mates  are  improved  significantly  after  smoothing.  In  ad¬ 
dition  it  is  found  that  the  effect  of  smoothing  can  be 
enhanced  if  the  computer  mapping  of  the  2D  model  field 
on  the  ID  vector  is  randomized:  Fig.  4a  bears  an  ap¬ 
parent  trace  of  columnwise  numbering  of  the  model 
field,  which  remains  visible  even  after  applying  the  al¬ 
gorithm,  generating  the  HM  columns.  As  a  consequence, 
error  distribution  in  Fig.  4a  contains  large  scales  in  cross¬ 
shore  direction,  which  tend  to  make  the  smoothing 
algorithm  less  effective  (Fig.  4c,  dashed  curves).  This 
drawback  can  be  easily  corrected  by  randomization  of 
the  above-mentioned  map  (the  randomized  HM  and  the 
RHM  method),  which  provides  an  error  pattern  similar 
to  Fig.  3a,  but  with  somewhat  lower  value  of  ( e )  clearly 
visible  in  Fig.  4c,  where  dashed  lines  show  evolution  of 
(e)  for  the  straight  HM  method  before  (upper  line)  and 
after  smoothing,  while  the  solid  black  lines  show  similar 
quantities  for  the  RHM  method. 

Comparison  with  the  MC  method  (gray  curves  in  Fig. 
4c)  demonstrate  a  noticeable  advantage  of  the  HM  tech¬ 
nique  (upper  curves),  which  remains  visible  at  higher  it¬ 
erations  n  >  100  even  after  smoothing  (lower  curves). 
This  advantage  increases  with  iterations  for  two  reasons: 
the  HM  method  converges  faster  than  fc~1/2  by  its  nature, 
whereas  the  efficiency  of  smoothing  (targeted  at  removing 
the  small-scale  error  constituents)  degrades  as  the  signal- 
to-noise  ratio  of  the  diagonal  estimates  increases  with  k. 

From  the  practical  point  of  view,  it  is  not  reasonable  to 
do  more  than  several  hundred  iterations,  as  (e)  drops  to 


FIG.  2.  Five  columns  of  the  BEC  operators  used  in  the  study:  (a) 
Gaussian-shaped  correlations  and  (b)  their  approximation  by  the 
inverse  of  the  second-order  polynomial  in  D.  White  circles  denote 
locations  of  the  diagonal  elements  of  the  corresponding  correlation 
matrices,  (c)  The  map  of  nonnormalized  diagonal  elements  of  Cp 
Depth  contours  are  in  m. 
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the  value  of  a  few  percent  (Fig.  3c),  which  is  much  smaller 
than  the  accuracy  in  the  determination  of  the  background 
error  variances.  We  may  therefore  conclude  that  it  is 
advantageous  to  use  the  RHM  technique  when  k  ~  100. 
In  this  case  (assuming  that  k  A/),  utilization  of  the  HM 
with  exact  dimension  is  not  necessary,  as  it  does  not  af¬ 
fect  the  convergence. 

d.  Asymptotic  expansion  method 

Since  the  principal  axes  of  the  diffusion  tensor  at  ev¬ 
ery  point  are  defined  by  construction,  computation  of 
the  zeroth-order  approximation  (6)  to  the  normalization 
factors  is  not  expensive.  Near  the  boundaries,  however, 
the  factors  described  by  (6)  have  to  be  adjusted  by  taking 
into  account  the  geometric  constraints  imposed  on  the 
diffusion.  In  this  study,  this  adjustment  was  computed 
for  points  located  closer  than  3A]  from  the  boundary.  It 
is  assumed  that  the  boundary  had  negligible  impact  on  the 
shape  of  the  diffused  8  function  (Fig.  5),  so  the  normali¬ 
zation  factor  near  the  boundary  was  computed  by  dividing 
the  reference  factor  (obtained  by  convolving  the  BEC 
operator  with  the  8  function  in  the  “open  sea”)  by  the 
respective  integral  over  the  land-free  subdomain  S'  shown 
in  Fig.  5: 

■,°(,)=2^[L.C(,’S,>8<X’y)H  '  (15) 


Integrals  in  the  rhs  of  (15)  have  to  be  taken  numeri¬ 
cally  for  all  the  near-boundary  points  x.  To  speed  up  the 
computations,  we  adopted  the  LH  assumption  near 
the  boundaries,  and  replaced  the  convolutions  in  (15)  by 
the  respective  analytical  functions  (4)-(5)  with  the  fixed 
value  i>(x).  It  is  necessary  to  note  that  the  assumption 
underlying  (15)  is  not  exact  for  the  zero  normal  gradient 
condition  in  use  (e.g.,  see  Mirouze  and  Weaver  2010). 
However,  the  errors,  caused  by  neglecting  distortions 
introduced  by  the  zero-flux  conditions  are  significantly 
smaller  (3%-7%,  see  Fig.  5)  than  the  accuracy  of  the  LH 
assumption  itself.  These  errors  could  be  removed,  for  ex¬ 
ample,  by  introducing  “transparent”  conditions  (Mirouze 
and  Weaver  2010). 

Figure  6  demonstrates  horizontal  distribution  of  the  er¬ 
ror  e(x)  obtained  by  approximating  the  diagonal  elements 
of  Cj  with  (6)  (i.e.,  the  zeroth-order  LH  method,  or  LH0) 


Fig.  3.  (a)  Error  distribution  after  60  iterations  of  the  MC  method, 
(b)  its  optimally  smoothed  version,  (c)  the  respective  dependences 
of  the  domain-averaged  error  <e)  and  on  the  iteration  number 
k .  The  dashed  line  is  the  approximation  of  by  the  k~m  law. 
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and  with  (13)  (i.e,,  the  first-order  LH  method,  or  LH1). 
Despite  an  apparent  violation  of  the  LH  assumption  in 
many  regions  (e.g.,  Aj  changes  from  205  to  the  back¬ 
ground  value  of  35  at  distances  L  ~~  5-65  <  A  i  across  the 
shelf  break),  the  mean  approximation  error  of  the  di¬ 
agonal  elements  appears  to  be  relatively  small  (19%)  for 
the  LHO  method,  with  most  of  the  maxima  confined  to 
the  regions  of  strong  inhomogeneity  (Fig.  6a).  The  next 
approximation  (Fig.  6b)  reduces  (e)  to  9%.  Numerical 
experiments  with  the  C2  model  have  shown  similar  re¬ 
sults  (16%  and  10%  errors). 

Another  series  of  experiments  are  performed  with  the 
varying  scaling  parameter  y  to  find  an  optimal  fit  to  d. 
Computations  were  made  for  0  <  y  <  L  The  best  result 
for  the  Gaussian  BEC  was  obtained  for  y2  =  0.30,  which 
is  fairly  consistent  with  the  value  (y2  =  0.33)  given  by 
(12).  In  the  case  of  C2  operator,  the  optimal  value  is  y2  = 
0.24,  still  in  a  reasonable  agreement  with  (12),  given  the 
strong  inhomogeneity  of  p  and  deviation  of  the  BEC 
operator  from  the  Gaussian  form.  A  somewhat  smaller 
value  of  y2(C2)  can  be  explained  by  the  sharper  shape  of 
the  respective  correlation  function  at  the  origin  (Fig.  1), 
which  renders  d°  to  be  less  dependent  on  the  inhomo¬ 
geneities  in  the  distribution  of  vy  and,  therefore,  requiring 
less  smoothing  in  the  next  approximation. 

e .  Numerical  efficiency 

Table  1  provides  an  overview  of  the  performance  for 
the  tested  methods.  For  comparison  purposes  we  show 
CPU  requirements  by  the  smoothed  MC  and  RHM 
methods  after  they  achieve  the  accuracies  of  the  LHO 
and  LH1  methods.  It  is  seen  that  both  MC  and  RHM 
methods  are  300-1000  times  more  computationally  ex¬ 
pensive  than- the  LH  technique.  In  fact,  for  the  2D  case 
considered,  the  computational  cost  of  the  stochastic 
methods  is  similar  to  that  of  the  three-dimensional 
variational  data  assimilation  (3DVAR)  analysis  itself. 

Comparison  of  the  spatial  distributions  of  the  ap¬ 
proximation  error  (e)(x)  favor  the  LH  methods  as  well: 
they  do  show  significantly  less  small-scale  variations  and 
may  have  a  potential  for  further  improvement. 

Comparing  Figs.  3b,  4b,  and  6b  shows  that,  in  contrast 
to  the  MC  and  RHM  methods,  LHO  errors  tend  to  increase 


Fig.  4.  (a)  Error  distribution  after  60  iterations  of  the  HM  method, 
(b)  As  in  (a),  but  for  the  smoothed  RHM  method,  (c)  Reduction  of 
the  domain-averaged  error  (e)  with  iterations  for  the  MC  (gray), 
RHM  (solid),  and  HM  (dashed)  methods.  The  bottom  graphs  are 
obtained  after  optimal  smoothing  of  the  diagonal  estimates.  Thin 
horizontal  lines  show  error  levels  provided  by  the  LHO  ((e)  =  0.17) 
and  LH1  methods  <e)  =  0.10. 
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in  the  regions  of  strong  inhomogeneity,  but  they  de¬ 
crease  substantially  after  smoothing  by  the  LH1  algo¬ 
rithm.  At  the  same  time,  the  LH1  errors  tend  to  have 
relatively  higher  values  near  the  boundaries;  the  effect  is 
less  visible  in  the  MC  and  RHM  patterns  (Figs.  3b  and 
4b).  This  feature  can  be  partly  attributed  to  certain  in¬ 
accuracy  in  the  algorithm  for  analytic  estimation  of  the 
near-boundary  elements  (Fig.  5c). 

Table  1  also  shows  that  LH  methods  outperform  both 
MC  and  HM  techniques.  Although  considerable  resources 
are  required  to  compute  near-boundary  integrals  for  the 
adjustment  factors  (15),  the  overall  CPU  savings  are 
quite  remarkable. 

/  LH  experiments  in  the  3D  setting 

To  check  the  performance  of  the  LHO  and  LH1  methods 
further,  a  larger  3D  domain  is  taken  from  the  NCOM 
setup  in  the  Okinawa  region  south  of  Japan  (Fig.  7),  with 
horizontal  resolution  of  10  km  and  45  vertical  levels. 
The  state  vector  dimension  M  (total  number  of  the  grid 
points)  in  this  setting  was  862  992. 

Because  of  the  large  My  it  is  computationally  un¬ 
feasible  to  directly  compute  all  the  diagonal  elements  of 
the  BEC  matrix.  Therefore,  accuracy  checks  are  per¬ 
formed  on  a  subset  of  10  000  points,  randomly  distrib¬ 
uted  over  the  domain  and  the  value  of  <e)  is  estimated  by 
averaging  over  these  points. 

The  diffusion  tensor  is  constructed  in  the  way  de¬ 
scribed  in  section  3a,  but  the  generating  field  v(x)  is 
taken  to  be  the  horizontal  velocity  field  from  an  NCOM 
run.  The  value  of  A3  (in  the  vertical  direction)  is  in¬ 
dependent  of  horizontal  coordinates,  but  varied  in  the 
vertical  as  3St,  where  Sz  is  the  vertical  grid  step.  Figure  7 
illustrates  spatial  variability  of  the  Cj  diagonal  elements 
at  z  =  20  m.  The  smallest  values  are  observed  in  the 
regions  of  Kuroshio  and  the  North  Equatorial  Current, 
where  the  largest  velocities  are  observed,  and  the 
fl  ~  \/detr  reaches  the  largest  values  [(6)].  To  better  test 
the  algorithm,  a  relatively  small  threshold  value  v  — 
0.02  m  is  prescribed,  so  that  diffusion  is  anisotropic 
in  more  than  90%  of  the  grid  points. 

Figure  8  demonstrates  the  accuracy  of  LHO  and  LH1 
methods  in  such  setting:  the  LHO  method  provides  an 


Fig.  5.  Adjustment  of  the  normalization  factors  near  the 
boundary,  (a)  Map  of  the  xth  column  of  C,  obtained  numerically  by 
integrating  (15)  over  the  real  domain  S'  (the  normalization  factor  is 
the  map  value  at  point  x).  (b)  Approximation  to  (a)  obtained  by 
integrating  Cj  in  the  infinite  domain  and  renormalizing  the  result  to 
have  the  same  integral  over  S'  as  in  (a),  (c)  The  difference  between 
(a)  and  its  approximation  (b). 


February  2012 


YAREMCHUK  AND  CARRIER 


645 


37.2 

36.0 

36.4 

36.0 

35  6 

~-123.2  -122.0  -122.4  -122.0  -121.6  -123.2  -122.8  -122.4  -122.0  -121.6 

Fig.  6.  Diagonal  approximation  errors  under  the  (a)  zeroth-order  and  (b)  first-order  LH  methods  for  the  Gaussian 

BEC  model. 


accuracy  of  9%,  which  is  further  improved  to  6%  by  the 
LH1  scheme.  The  major  improvement  occurs  in  the  re¬ 
gions  where  points  with  highly  anisotropic  v  neighbor 
isotropic  points  and  reduce  the  diagonal  elements  in 
the  latter.  The  effect  is  reflected  by  the  negative  bias 
of  the  scatterplot  at  high  values  of  d°,  which  reach 
their  maximum  of  0.0237  in  the  points  with  isotropic  v 
(Fig.  8a). 

Figure  8c  shows  the  dependence  of  approximation  error 
e  on  the  value  of  y3  for  both  correlation  models.  The  best 
approximation  is  obtained  at  y3  =  0.26,  a  value  somewhat 
lower  than  suggested  by  the  heuristic  formula  (y3  - 
5/18  =  0.28).  Similarly  to  the  2D  case,  the  optimal  value 
of  y3(C2)  =  0.21  is  less  than  y3(C1),  in  agreement  with 
more  rapid  off-diagonal  decay  of  the  C2  matrix  elements. 

In  general,  it  appears  that  the  relationship  (12)  pro¬ 
vides  a  reasonable  guidance  to  the  estimation  of  the 
smoothing  parameter  in  the  LH1  method.  For  the  C, 
model,  the  operator  acting  on  d°  can  be  implemented  by 
either  reducing  y'1  times  the  number  of  time  steps  in  the 
integration  of  the  diffusion  equation,  or  by  y~1/2-fold 
reduction  of  the  decorrelation  radius.  For  the  C2  model 
only  the  second  option  is  applicable:  it  would  also  reduce 
the  number  of  iterations  required  for  computing  the  ac¬ 
tion  of  the  BEC  operator. 


Table  1.  Relative  CPU  times  required  by  the  MC  and  RHM 
methods  to  achieve  the  accuracies  <s)  of  the  LH0  and  LH  1  methods 
(shown  in  parentheses). 


MC/LH0 

MC/LH1 

RHM/LH0 

RHM/LH1 

C,  755  (0.19) 

1205  (0.09) 

680  (0.19) 

520  (0.09) 

C2  780  (0.17) 

490  (0.10) 

850  (0.17) 

330  (0.10) 

4.  Summary  and  discussion 

In  this  study  we  examined  the  computational  effi¬ 
ciency  of  several  techniques  used  for  estimating  the  di¬ 
agonal  elements  of  the  two  types  BEC  operators:  with 
the  Gaussian-shaped  kernel  C{  and  with  the  kernel  gen¬ 
erated  by  the  second-order  polynomial  approximation  to 
Cr  The  considered  techniques  include  the  “stochastic” 
MC  and  HM  methods,  which  retrieve  diag(C)  only 
from  its  action  on  a  vector,  and  the  “deterministic” 
scheme  based  on  the  analytic  diagonal  expansion  under 
the  assumption  of  local  homogeneity  of  the  diffusion 
tensor.  The  deterministic  scheme  was  tested  in  two 


Fig.  7.  Diagonal  elements  of  Cj  in  the  Okinawa  region  at  z  =  20  m. 
The  actual  values  are  multiplied  by  104. 
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regimes:  the  zeroth  (LHO)  and  the  first-order  (LH1) 
approximations. 

Numerical  experiments  conducted  with  realistic  dif¬ 
fusion  tensor  models  have  shown  that  (i)  HM  technique 
proves  to  be  superior  in  efficiency  compared  to  the  MC 
technique  when  accuracies  of  less  than  10%  (fc  >  100) 
are  required;  (ii)  both  stochastic  methods  require  300- 
1000  times  more  CPU  time  to  achieve  the  accuracy,  com¬ 
patible  with  the  most  efficient  LH1  method;  (iii)  with  the 
Gaussian  model,  the  LH1  method  demonstrates  the  best 
performance  with  the  value  of  the  smoothing  parameter 
y  compatible  with  the  one  given  by  the  relationship  (12) 
derived  from  the  asymptotic  approximation  of  the 
Gaussian  kernel  diagonal. 

In  deriving  the  ansatz  (13)  for  the  LH1  model  we  fol¬ 
lowed  the  approach  of  Purser  et  al.  (2003),  who  pro¬ 
posed  to  smooth  the  zeroth-order  diagonal  by  the  square 
root  of  the  BEC  operator  in  a  one-dimensional  case.  Using 
the  asymptotic  technique  for  the  heat  kernel  expansion,  we 
obtained  a  formula  for  higher  dimensions,  and  tested  its 
validity  by  numerical  experimentation. 

It  should  be  noted  that  the  formal  asymptotic  expan¬ 
sion  (7)  is  local  by  nature  and  tends  to  diverge  in  practical 
applications,  where  spatial  variations  of  the  diffusion 
tensor  may  occur  at  distances  L  comparable  with  the 
typical  decorrelation  scale  A.  To  effectively  immunize 
the  expansion  from  the  ill  effects  of  the  abrupt  changes 
in  v,  we  utilized  a  nonlocal  empirical  modification,  still 
fully  consistent  with  the  original  expansion  in  the  limit 
A/L  -►0,  but  sufficiently  robust  with  respect  to  the  errors 
related  to  the  high-order  derivatives  of  v.  A  similar 
technique  was  developed  by  Purser  (2008a, b),  who  used 
empirical  saturation  functions  to  stabilize  the  higher- 
order  approximations  of  the  heat  kernel  diagonal. 

In  general,  results  of  our  experiments  show  high  com¬ 
putational  efficiency  of  the  LH1  scheme,  whose  total  CPU 
requirements  is  just  a  fraction  of  the  CPU  lime  required 
by  the  convolution  with  BEC  operator — a  negligible 
amount  compared  to  the  cost  of  3D  VAR  analysis. 

A  separate  question,  requiring  further  investigation, 
is  the  accurate  treatment  of  the  boundary  conditions.  In 
the  present  study  we  assumed  that  boundaries  affect  only 
the  magnitude  of  the  corresponding  columns  of  C,  but 
not  their  structure.  This  approximation  is  only  partly 


Fig.  8.  Scatterplots  of  the  true  diagonal  elements  of  C1  (vertical 
axis)  vs  their  approximations  by  (a)  LHO  and  (b)  LH1  algorithms. 
The  actual  values  are  multiplied  by  103.  Near-boundary  points  are 
excluded,  (c)  Diagonal  approximation  errors  as  a  function  of  y  for 
the  (black)  and  C2  (gray)  BEC  models.  The  dashed  line  shows 
the  value  of  y3  given  by  (12). 
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consistent  with  the  zero  normal  flux  conditions  for  D,  but 
can  be  avoided  if  one  uses  “transparent”  boundary  con¬ 
ditions  (e.g.,  Mirouze  and  Weaver  2010),  which  do  not 
require  computation  of  the  adjustment  factors  (15). 

On  the  other  hand,  it  might  be  beneficial  to  keep 
physical  boundary  conditions  in  the  formulation  of  D,  as 
they  are  likely  to  bring  more  realism  to  the  dynamics  of 
the  BE  field.  In  the  considered  diffusion  tensor  model, 
anisotropic  BE  propagation  is  superimposed  on  the 
small-scale  isotropic  BE  diffusion,  which  takes  place  at 
scales  that  are  not  well  resolved  by  the  grid  (less  than 
35).  This  may  give  some  grounds  to  employ  isotropic 
homogeneous  model  at  distances  d  <33  from  the  bound¬ 
ary,  because  the  deterministic  part  of  the  BE  transport 
associated  with  the  boundary  effects  is  not  well  resolved 
at  these  scales  anyway.  Such  simplification  greatly  re¬ 
duces  the  computational  cost  of  the  LH  algorithm  with 
zero-flux  conditions  because  the  region  screened  by  the 
boundary  is  readily  available  from  the  land  mask  and 
needs  no  computation/rescaling  at  every  point.  In  con¬ 
trast  to  the  2D  problem  considered,  such  computations 
may  be  more  costly  in  the  3D  problems,  as  the  number 
of  adjustment  factors  n  ~  may  be  quite  large  even 
with  a  moderate  grid  size  N  ~  106-107.  However,  our 
preliminary  experiments  with  inhomogeneous  BEC 
models  with  larger  N  have  shown  that  estimation  of  the 
diagonal  with  the  LH1  scheme  still  remains  two  orders  in 
magnitude  less  expensive  than  the  3DVAR  analysis. 

Results  of  this  study  indicate  that  LHI  approximations 
to  the  BEC  diagonal  may  serve  as  an  efficient  tool  for 
renormalization  of  the  correlation  operators  in  variational 
data  assimilation,  as  they  are  capable  of  providing  3%- 
10%  accuracy  in  realistically  inhomogeneous  BEC  models. 
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APPENDIX  A 
Hadamard  Matrices 

By  definition,  a  Hadamard  matrix  (HM)  is  a  square 
matrix  whose  entries  are  either  1  or  -1  and  whose  col¬ 
umns  are  mutually  orthogonal.  The  simplest  way  to  con¬ 
struct  HMs  is  the  recursive  Sylvester  algorithm,  which 
is  based  on  the  obvious  property:  if  is  an  M  X  M 
Hadamard  matrix,  then 


is  also  HM.  Starting  from  H2  =  [1  1;  1  —1],  the  HMs  with 
order  M  =  2",  n  =  1,2...  can  be  easily  constructed.  The 


HMs  with  M  =  12, 20  were  constructed  “manually”  more 
than  a  century  ago.  A  more  general  HM  construction 
algorithm,  which  employs  the  Galois  fields  theory  was 
found  in  1933.  However,  it  is  still  unknown  if  HMs  exist 
for  all  M  =  An  where  n  is  a  positive  integer. 

Hadamard  matrices  are  widely  used  in  many  branches 
of  applied  mathematics  and  statistics  (http://en.wikipedia. 
org/wiki/Hadamard_matrix).  In  the  present  study  we 
used  the  MatLab  code  that  handles  only  the  cases  when 
M/12  or  Ml 20  is  a  power  of  2.  Despite  this  restriction,  the 
available  values  of  M  were  sufficient  for  our  purposes. 

APPENDIX  B 

Correlation  Modeling  with  Diffusion  Operator 

Consider  a  family  of  linear  operators  B  in  the  space  7i 
of  sufficiently  smooth  square-integrable  functions  /(x), 
x  €  Rn: 

(I  -  D/2m)mB(x,x')  =  S(x,x').  (Bl) 

Assume  that  v  =  const  in  the  definition  of  the  diffusion 
operator  (1)  and  perform  the  coordinate  transformation 
x  =  p~mx.  In  the  new  coordinate  system  the  diffusion 
tensor  v  is  the  unity  matrix,  D  takes  the  form  of  the 
Laplacian  operator  A,  whereas  the  3  function  in  the  rhs  of 
(Bl)  is  rescaled  as  a  result  of  the  contraction  fl  =  \Zdety 
of  the  volume  element:  6(x,x')  fTl5(x,x').  Equation 

(Bl)  takes  the  following  form: 

(1  -  A/2/n)mB(x,x')  =  (B2) 

After  changing  the  coordinates  in  H  with  the  Fourier 
transform,  both  A  and  the  operator  B  are  diagonalized. 
The  diagonal  elements  are  given  by 

B(k)  =  ^(1  +  k2/2m)~m.  (B3) 

It  is  easy  to  note  that  by  virtue  of  the  textbook  formula 
exp*  =  lim^^l  +j t/w)m,  the  Fourier  representation 
(B3)  of  B  converges  at  large  m  to 

B1(k)=iexp(-k2/2),  (B4) 

the  well-known  expression  for  the  Green  function  of  the 
diffusion  equation.  Since  the  diagonals  given  by  (B3)- 
(B4)  are  strictly  positive,  the  respective  operators  in  H 
are  positive  definite,  and  can  be  interpreted  as  correla¬ 
tion  operators.  This  property  has  been  widely  used  for 
the  background  error  covariance  modeling  in  practical 
data  assimilation. 
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Matrix  dements  of  the  operators  =  exp(D/2)  and 
B2  =  (I  —  D/2 m)~m  can  be  found  via  the  inverse  Fourier 
transformations  of  (B3)-(B4): 

Bi0)  =  |2^}r„  B,(k)exp(ikr)<fk 


B2(r)  =  (2^  jR„  B(k)  expOkr)  dk 

1  \\f2mt\s  Ks(\\/2mr\) 

(2ir)nf2n  2-ir  (5)  ’ 

where  s  =  m  -  nl2  and  r2  =  (x{  —  Xj)1'^  -  x2)  1S 
square  of  the  distance  between  the  correlated  points  Xj 
and  itj  in  the  transformed  coordinate  system.  Making  a 
substitution  r2  =  (Xj  —  x2)Tv“1(x1  —  x2)  in  the  rhs  of 
(B5)-(B6),  we  obtain  the  matrix  elements  (4)-(5)  in  the 
original  coordinate  system. 

From  the  viewpoint  of  numerical  applications  it  is  in¬ 
structive  to  connect  the  discretized  operator  equation  (Bl) 
with  various  time  integration  techniques  of  the  dis¬ 
cretized  diffusion  equation.  The  numerical  approxima¬ 
tion  of  B  in  (Bl)  is  never  calculated  in  practice  because 
of  the  immense  cost  of  such  a  computation.  Instead,  the 
result  /™(x)  of  action  by  B  on  a  (discrete)  model  state 
vector  /°(x)  is  calculated  by  solving  the  system  of 
equations 

(I  -  D/2m)mfm  =  f°,  (B7) 

where  D  denotes  the  discretized  diffusion  operator.  If 
we  now  assume  that/°(jc)  represents  the  “initial  state” 
and  prescribe  the  “time  step”  5/  such  that  the  “virtual 
ntegration  time”  is  m8t  —  1,  the  action  of  the  correlation 
operator  (Bl)  can  be  identified  as  a  result  of  discrete-time 
integration  of  the  diffusion  equation  BJ  =  D/2/  with  the 
implicit  scheme: 

/'■  -  rX  =  \stbf,  (B8) 

starting  from  the  initial  state/0. 

Similarly,  the  action  of  exp(D/2)  is  never  computed  by 
convolving  a  state  vector  f°  with  the  discretized  kernel 
(B5),  but  rather  by  the  discrete-time  integration  of  the 
diffusion  equation  with  the  explicit  numerical  scheme: 

P  -  1  =  \stbr\  i  =  1,...  ,m.  (B9) 


In  contrast  to  the  implicit  scheme  (B8),  5/  in  this  numerical 
method  is  limited  from  above  by  the  stability  condition 
requiring  that  eigenvalues  of  the  operator  I  +  5/6/2  must 
be  less  than  1  in  the  absolute  value.  As  a  consequence, 
the  minimum  number  of  time  steps  m  =  1/5/  may  be  quite 
large,  making  the  numerical  implementation  of  exp(D/2) 
computationally  impractical. 
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